Monte Carlo simulation of uncoupled continuous-time random walks yielding a 
stochastic solution of the space-time fractional diffusion equation 
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We present a numerical method for the Monte Carlo simulation of uncoupled continuous-time 
random walks with a Levy a-stable distribution of jumps in space and a Mittag-Leffler distribution 
of waiting times, and apply it to the stochastic solution of the Cauchy problem for a partial differ- 
ential equation with fractional derivatives both in space and in time. The one-parameter Mittag- 
Leffler function is the natural survival probability leading to time-fractional diffusion equations. 
Transformation methods for Mittag-Leffler random variables were found later than the well-known 
transformation method by Chambers, Mallows, and Stuck for Levy a-stable random variables and 
so far have not received as much attention; nor have they been used together with the latter in spite 
of their mathematical relationship due to the geometric stability of the Mittag-Leffler distribution. 
Combining the two methods, we obtain an accurate approximation of space- and time-fractional 
diffusion processes almost as easy and fast to compute as for standard diffusion processes. 

PACS numbers: 02.50.Ng, 02.70.Tt, 02.70.Uu, 05.70.Ln 

I. INTRODUCTION II. THEORY 



Continuous-time random walks (CTRWs) and frac- 
tional diffusion equations (FDEs), or fractional Fokker- 
Planck equations, have received increasing attention. 
Metzler and Klafter reviewed analytical and numerical 
methods to solve fractional equations of diffusive type 
Q. In Refs. 0, i, i, H, i, 0, i, i, 0, applications and 
enhancements of these techniques were presented. The 
relevance of fractional calculus in the phenomenologi- 
cal description of anomalous diffusion has been discussed 
within applications of statistical mechanics in physics, 
chemistry and biol2gyilMl [H, Q E [M [13 as well 
as finance [H, [l^, [lo, l2lU22j|: even human travel and the 
spreading of epidemics were modeled with fractional dif- 
fusion [231 . A direct Monte Carlo approach to fractional 
Fokker-Planck dynamics through the underlying CTRW 
requires random numbers drawn from the Mittag-Leffler 
distribution. Since sampling the latter was considered 
troublesome, different schemes to avoid it were proposed. 
One possibility consists in replacing it with the Pareto 
distribution — i.e., its asymptotic power-law approxima- 
tion for t — > oo [24|; however, as the authors point out, 
this is limited to long times and an index (3 not close to 
1. A more general alternative is based on subordination 
[25I . [26l . Here we present a straightforward Monte 

Carlo method for the efficient simulation of uncoupled 
CTRWs using an inversion formula for the Mittag-Leffler 
distribution and apply it to compute approximate solu- 
tions of the Cauchy problem for a generalized diffusion 
equation that has fractional space and time derivatives. 
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A. Continuous-time random walks 

A CTRW 28] is a pure jump process; it consists of a se- 
quence of independent identically distributed (i.i.d.) ran- 
dom jumps (events) separated by i.i.d. random waiting 
times Ti, 

n 

tn = ^n, neR+, (1) 

SO that the position at time t € [tmtn+i) is given by 

n 

x{t) = Y.^,, eR. (2) 

1=1 

A realization of the process is a piecewise constant func- 
tion resulting from a sequence of up or down steps with 
different height and depth; see Fig.[TJ Jumps are assumed 
to happen instantaneously or at least in negligible time. 
In general, jumps and waiting times depend on each other 
and they can be described by a joint probability density 
(y9(^,r). The latter appears in the integral equation giv- 
ing the probability density p(x, t) for the process being 
in position x at time i, conditioned on the fact that it 
was in position a; = at time t — Q: 

/ + 00 I't 
/ dT^{tT)p{x-^,t-T). 
-cc Jo 

(3) 

Here the initial condition a;(0) = is contained implicitly 
in the first term S{x)'^{t), where we find the complemen- 
tary cumulative distribution function (survival function) 

/ + 00 rt 
/ dr^(^,T). (4) 
oo Jo 

Recently, one of the authors of this paper presented an 
analytical solution of the integral equation in the uncou- 
pled case — i.e., when 'f{S,,T) = X{^)iP{t), where A(^) is 
the jump marginal density and iIj(t) is the waiting time 
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FIG. 1: (Color online) Sample paths of CTRWs with scale 
parameters 7t = 0.001, 7^, = 7f''", and different choices of a 
and 13. With smaller a the jumps become larger; with smaller 
P the waiting times become longer. 



B. Fractional diffusion equation 



The well-known standard diffusion equation 



d d'^ 
Q-^u{x,t) = D—u{x,t), 

u{x,0+) = d{x), xeR, te 



(5) 



can be generalized to the space-time fractional diffusion 
equation 

g^u{x,t) = D^^u{x,t) (6) 

u{x,0+) = S{x), xeR, teR+, 

where, for < a < 2, 9"/9|a;|" denotes the symmetric 
Riesz-Feller operator of symbol — and, for < (3 < 1, 
/dt^ is the Caputo derivative [IM Hi] . Without loss 
of generality, we assume D — 1] a. different value would 
just mean a scale transformation of space and/or time 
units. u{x,t) > is the Green function of the FDE, 

uix, t) = r^'°' W{x/tP/''; a, (3), (7) 

with the scaling function 

W{ta,P)^T-'[Ef,{-\nn]{0- (8) 
Ep{z) is the one-parameter Mittag-Leffler function [s^ . 



n=0 



r(/3n + 1) ' 



z e C, 



(9) 



with 



(t), teR+. (10) 



_l + s/3_ 

and C denote the Fourier and Laplace transforms: 



-l-oo 



fi^)=^Afi^)m= f{x)e^^^dx, (11) 

^ — oo 

f(t)e-'Ut. sec. (12) 



For t gM. and (3=1, the Mittag-Leffler function with ar- 
gument —t^ reduces to a standard exponential decay e~*; 
when < (3 < 1, the Mittag-Leffler function is approxi- 
mated for small values of t by a stretched exponential de- 
cay (WeibuU function) exp{—t^/a), where a = r(/3 -I- 1), 
and for large values of t by a power law bt~^, where 
b = r(^)sin(/37r)/7r; see Fig. EI The Mittag-LefHer dis- 
tribution is an important example of fat-tailed waiting 
times; it arises as the natural survival probability leading 
to time-fractional diffusion equations. There is increasing 
evidence for physical phenomena [H, [s^, [IHl and human 
activities [S^, [33, [111 that do not follow either exponential 
or, equivalently, Poissonian statistics. 
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FIG. 2: (Color online) The Mittag-Leffler complementary cu- 
mulative distribution function sampled from Eq. (|20p (circles) 
and computed analytically (solid line) [39l ]. as well as its ap- 
proximations for f — > (WeibuU function, long dashes) and 
t ^ oo (power law, short dashes). 

Equations ^ and ^ can be obtained by Fourier- 
Laplace transformation of the FDE, recalling the defi- 
nition of the fractional derivatives used in Eq. 

The space- fractional derivative of order a G (0,2] is 
defined according to Riesz [40l ]: 



fix)^T-' -iK^fiK) (x) 



(13) 



For a = 2 this reduces to the usual second order deriva- 
tive. For a < 2 the following equation holds: 



d'^f(x) r(a + l) an 

„ , — = sm — 

d|a;|" TT 2 



' f{x+^)^2f{x) + f{x~0 



(14) 

The time-fractional derivative of order /3 E (0, 1] is 
defined according to Caputo [4l|, \^ : 



dtp 



/(0+) it). (15) 



For (3 = 1 this reduces to the usual first order derivative. 
For /? < 1 the following equation holds: 



dt0 



1 



r{i-(3) 



dt 



fir) 



it 



■dT 



/(0+) 



, (16) 



where /(O"*") is the initial condition. For a = 2 and (3=1, 
the standard diffusion equation, Eq. ([5]), is recovered. 
It is inevitable to solve numerically a FDE in the most 
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equation, which may include space- and time-dependent 
diffusion and drift terms. Possible approaches are the 
direct calculation of the integrals in Eq s. (fT4|) and (fT6|) 
[4^ , finite-difference methods [H, |4y| , and stochastic 
methods ^5,, J, [13, [H, [13 . All of them are complicated, 
the latter ones mainly because of the supposedly cum- 
bersome generation of Mittag-Leffler random numbers. 
While this problem has been often worked around in the 
past, we show how to overcome it, obtaining a fast and 
accurate method for the Monte Carlo solution of FDEs 
via uncoupled CTRWs. As a benchmark, we focus our 
attention on the Cauchy problem defined in Eq. ([5]), for 
which an analytical solution given by Eqs. (O and ([8]) is 
available. 



C. Link between continuous-time random walks 
and the fractional diffusion equation 

The link between CTRWs and time-fractional diffu- 
sion was discussed rigorously in Ref. fi?! in terms of the 
generalized Mittag-Leffler function Ep_p{—T^). 

In order to approximate the Green function in Eq. ([7]), 
it is sufficient to simulate CTRWs whose jumps are dis- 
tributed according to the symmetric Levy a-stable prob- 
ability density (which reduces to a Gaussian for a = 2) 



i„(0=^-Mexp (-17,^1")] (0 



(17) 



and whose waiting times have the probability density 



V'/3(t) 



dr 



(18) 



where Ep[z) is the one-parameter Mittag-Leffler function 
given by Eq. ([9]). Then a weak-limit approximation of the 
Green function is obtained by rescaling waiting times by 
a constant 74 and jumps by a constant = 7f^": let- 
ting 74 (and as a consequence 7^;) vanish, and plotting 
the histogram for the probability density p.^^ ^j^{x,t; a, /?) 
of finding position x at time t for the rescaled process. 
This probability density weakly converges to the Green 
function u(x,t] a, (3). Weak convergence means that for 
a; = a singularity is always present in p~f^,^^{x,t; a,/9) 
at a; = for any finite value of jt and ■ This singular- 
ity is the term (5(a;)*(t) in Eq. §1 with «'(t) = Ep{~tf^). 
In the case a = 2 and (3=1 the CTRWs are normal 
compound Poisson processes (NCPPs) and, in the diffu- 
sive limit, one recovers the Green function for the stan- 
dard diffusion equation, Eq. ([5]) — i.e., the Wiener pro- 
cess. This procedure is justified in Refs.[l and] 2 9.. In the 
latter reference, one can also find a theoretical justifica- 
tion for the Monte Carlo procedure where waiting times 
are generated according to a power-law distribution; a 
more complete treatment has been given in Ref. i25l . 



III. TRANSFORMATION FORMULAS FOR 
NON UNIFORM RANDOM NUMBERS 

The usual methods for generating random numbers 
with a specific probability density are transformation, 
also called inversion because it requires the inverse cu- 
mulative distribution function [48j . and von Neumann 
rejection [4^. While the latter is more general, the for- 



Symmetric Levy a-stable probability 
distribution 



The symmetric Levy a-stable probability density 
LaiO foi' the jumps, Eq. ([T7)) . can be calculated by series 
expansion, which we do not report here, by direct inte- 
gration [50I . [51I or by numerical Fourier transform . 
These methods produce a pointwise representation of the 
density on a finite interval that can be used for rejection, 
most efficiently with a lookup table and interpolation. 
More convenient is the following transformation method 
by Chambers, Mallows, and Stuck [ssj : 



Co 



— log U COS 4> 

cos((l — a)0) 



sai{aq 



(19) 



where (p = 7r(w — 1/2), u,v € (0, 1) are independent uni- 
form random numbers, is the scale parameter, and 
is a symmetric Levy a-stable random number. For 
a = 2, Eq. (fTOj) reduces to ^2 = "^^xy/— log u sin 0, i.e. 
the Box-MuUer method for Gaussian deviates. The other 
two notable limit cases are the Cauchy distribution, with 
a = 1 and ^1 = jx tan(/), and the Levy distribution, with 
a = 1/2 and ^1/2 = —7a; tan 0/(2 log u cos ^). 



B. One-parameter Mittag-Leffler probability 
distribution 



The probability density tppiT) for the waiting times, 
Eq. (jlSp . can be computed as a power series from the 
definition of the one-parameter Mittag-Leffler function, 
Eq. ([9]), leading to a pointwise representation on a fi- 
nite interval; random numbers can then be produced 
by rejection, again with a lookup table and interpola- 
tion. Though CTRW sample paths with a Mittag-Leffler 
waiting time distribution have appeared in the literature 
[25, 26; 27, 54], so far it has not been recognized in this 
context that inversion formulas analogous to Eq. are 
available [H, [H, [13, [H, [I§, [S^, [Ml, [12 ■ The most conve- 
nient expression is due to Kozubowski and Rachev |5l| : 



f sm{(3n) ^ 
Tf3 = -7t logu ; — -5 — r - cos(/37r) 
\tan(/57rw) 



(20) 



where u,v € (0, 1) are independent uniform random num- 
bers, 7t is the scale parameter, and is a Mittag- 
Leffler random number. For /3 = 1, Eq. ([20]) reduces 
to the inversion formula for the exponential distribution: 
Ti = — 7t logM. Equation (|20p and equivalent forms stem 
from mixture representations of a Mittag-Leffler random 
variable through an exponential and a stable random 
variable. The oldest representation is [sl, [6l| 



1//3 



(21) 



where is a skew Levy a-stable random number inde- 
pendent of Ti, with index a = (3, skewness parameter 1, 
and scale factor jx = 1/8. A more recent representation 

is MM 



Tf) = Tl 4 



±1/13 



(22) 



where ^1+ is a positive random number distributed ac- 
cording to a Cauchy distribution Li+(^) with scale pa- 
rameter jx = sin(/37r), location parameter S = — cosiPn), 
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The connection of Mittag-LefHer to stable random vari- 
ables can be obtained in the framework of the theory of 
geometric stable distributions. A random variable ^ is 
stable if and only if, for all n e N i.i.d. copies of it, 
^1, . . . , there exist constants a„ G K+ and fo„ G K such 
that the scaled and shifted sum a„ (^i + • • • + + 6„ has 
the same distribution as ^. A Mittag-Leffler random vari- 
able is not stable, but it is geometric stable [i^; i.e., it is 
the weak limit for p of the appropriately scaled and 
shifted geometric random sum a{p) [ti -|- • • • -I- Tj,(p)] -I- b{p) 
of suitable i.i.d. random variables Tj, where iy{p) is a geo- 
metric random variable indepedent of each r,; , with mean 
p € (0, 1), and a geometric probability distribution 

P{i,{p) = n) ^ p{l - p)"-\ neN. (23) 

A random variable is geometric stable if and only if its 
characteristic function V-'('*) is related to the character- 
istic function A(k) of a stable random variable by the 
equation [6^] 

= , (24) 
1 - logA(K) 

With this one-to-one correspondence, a parametrization 
of a geometric stable probability density can be es- 
tablished from a parametrization of the corresponding 
stable probability density X{x). Geometric random sums 
of symmetric Ti yield the class of Linnik distributions (a 
generalization of the Laplace distribution ^e"'*'), while 
positive Ti yield the class of Mittag-Leffler distributions 
(as already seen, a generalization of the exponential dis- 
tribution e~*, t > 0). In particular, the Mittag-LefHer 
distribution can be written as a mixture of exponential 
distributions (iil,^]: 



Ef}{-t'^)= / exp(-^i).g(/i) d/i, 
Jo 

with a weight 

1 sin(/37r) 



9{^A = 



TT fi^+l^^ + 2 cos(/37r)/i + 

given by g{ii)dfi — Li+{iil^)dii^ , where Li+(^) is the 
probability density of in Eq. (|22p introduced be- 
fore. Equations (|25p and ([^5]) express Eq. in terms 
of density functions. The inverse cumulative distribu- 
tion of ii+(^) yields the transformation formula for 
appearing as the argument of the power function in 
Eq. ([20|) [ssl . Isgj ]. Alternatively, the inversion formula 
^1 = 7a;tan(/() -|- 5 for ii(C), see Eq. (ITOl) . can be substi- 
tuted into Eq. ([22]) . provided negative values of are 
discarded. 

An older equivalent form of Eq. ([20|l was obtained sub- 
stituting an inversion formula for [Q§\ into Eq. (PT|) 
jssl . [6l| . A similar result can be reached using a general 
transformation formula for skew Levy a-stable random 
numbers [s^l, of which Eq. (fT9|) is a special case with 
skewness parameter 0. Both ways require three indepen- 
dent uniform random numbers and more transcendent 
functions than Eq. ([20| . making the latter slightly more 
appealing from a numerical point of view. 



(25) 



(26) 



IV. NUMERICAL RESULTS 



Examples of CTRWs generated according to the de- 
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TABLE L Average number n of jumps per run and total CPU 
time tcpu in seconds for 10^ runs with t G [0, 2] on a 2.2 GHz 
AMD Athlon 64 X2 Dual-Core with Fedora Core 4 Linux, us- 
ing the rani uniform random number generator [67l | and the 
Intel C-|— I- compiler version 9.1 with the -03 -static optimiza- 
tion options. 



shown in Fig. [TJ The complementary cumulative distri- 
bution function (survival function) of random numbers 
obtained through Eq. (|20p is checked against its analytic 
value |39] and its approximations for t Q and t — *■ c» 
in Fig. [51 where a log-log scale and logarithmic binning 
[6^ is used. Timings are reported in Table |T] and Ref. IgS 
The advantage of Eq. ([20]) is that Mittag-Leffler devi- 
ates are generated with a simple and elegant procedure 
and no accuracy losses due to truncation of the power se- 
ries in Eq. ([9]) or truncation of the density function to a fi- 
nite interval as necessary in the rejection method. The ef- 
fects of the truncation of the jump density in Levy flights 
are analyzed in Ref . [H, whereas no study is available for 
truncation effects on Mittag-Leffler deviates. Together 
with Eq. (Uni), a scheme is obtained that yields sample 
paths for a CTRW with a Levy jump marginal density 
and a Mittag-Leffler waiting time marginal density at a 
speed comparable to that of a NCPP: Though each point 
for a generic CTRW takes about 3.6 times more than for 
a NCPP, fewer points are necessary (see n in Table |T| be- 
cause the waiting times are longer. The latter reference 
reports also that if Levy and Mittag-Leffler random num- 
bers are produced by rejection, computing the values of 
the probability density functions simple-mindedly with a 
series expansion every time they are needed, rather than 
just once at the beginning to set up a lookup table, for 
Levy deviates the procedure takes 400 times longer than 
with Eq. and for Mittag-Leffler deviates it takes 5000 
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FIG. 3: (Color online) Decay of the probability density 
Pj^,jt{x,t; a, f3) with a — 1.7, (3 — 0.8, jt ~ 0.1, and 
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FIG. 4: (Color online) Convergence of t^^°'p-/^,-/t{x,t; a, (3) 
to the scaling function W{x/t'^^°'; a,/3), Eq. (O, at t = 2 for 
selected values of a and /3. The curves are shown in a time- 
independent way as scaling plots and appear in the same or- 
der from bottom to top as reported in the legend — i.e., with 
decreasing 74. The curve with the smallest jt is almost in- 
distinguishable from its theoretical limit W (solid black line) . 
However, in spite of the impression that may arise from the 
few terms and the ranges chosen here, in general the function 
sequences are not monotonic. The scale parameters 7^, and 
7t tend to as 7" = 7f . The central peak decreases when 
the ratio t/^t becomes larger, as is evident in Fig. [3] 



times longer than with Eq. (|^D|) . Because of the slow con- 
vergence of the power series in Eq. up to 200 terms 
are necessary to achieve an acceptable accuracy, and each 
term is computationally expensive because of the F func- 
tion. Of course these are extreme figures on the other 
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FIG. 5: (Color online) Convergence of max3;=^o |P7x,7t (2;, 
a, 13) — u{x,t; a,/3)| for selected values of a and (3 when 
yx, 7t ^ with 7" = 7f . 



latter can be; there are smarter ways to compute both 
the Levy and Mittag-LefHer [S^H^I probability densities. 

Using many CTRW reahzations, histograms can be 
built that give the evolution of p{x, t) with initial con- 
dition p{x,0) = S{x), as displayed in Fig. [3l Accord- 
ing to Eq. ([3]), the initial condition evolves as S{x)'^{t); 
i.e., it is visible as a spike at a; = that decays as t 
evolves. The mass of the spike is 4'(i) = Efj {—{t/jt)^)- 
In Fig. [3] this feature appears as a crest. Figure |4] shows 
how histograms built with CTRWs converge to the Green 
function, Eq. ([7]), of the FDE for decreasing values of 

the scale parameters 74 and — lt^°'- To evaluate 
the scaling function in Eq. ([8]) needed for Eq. jTl), we 
used standard algorithms for Epi^—t^) [s^, [39|, in- 
cluding the fast Fourier transform. In Fig. [5] we plot 
maxa;^o b7x,7t(^i^; Q^,/3) — u{x,t; a, l3)\ as a function of 
vanishing 74 with "f^ = "ft^" ■ A rigorous analysis of con- 
vergence bounds is beyond the scope of this paper. 



V. CONCLUSIONS 

The use of Mittag-Leffler random numbers generated 
according to Eq. (|20[) in combination with Levy random 
numbers generated according to Eq. (|19p is very useful 
in the Monte Carlo simulation of uncoupled continuous- 
time random walks. In the hydrodynamic limit, ap- 
propriately rescaled uncoupled continuous-time random 
walks with a one-parameter Mittag-Leffler distribution of 
waiting times and a symmetric Levy a-stable distribution 
of jumps in space yield the Green function of the Cauchy 
problem for a space-time fractional diffusion equation; 
we verified this for Eq. which has an analytical so- 
lution, Eq. ([7]), as a benchmark for more difficult cases 
where the diffusion and drift terms depend on space and 
time. We have shown that the computational effort for 
a fractional diffusion process is almost as small as for a 
standard diffusion process. It is true that in the same 
fluid limit the Green function can be obtained too by 
Monte Carlo sampling of just the asymptotic power-law 
tail approximations of the Levy and Mittag-Leffler prob- 
ability distributions, at least when the indices a and /3 
are not close to 2 and 1, respectively. However, the neat 
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numerically so convenient that there is no good reason for 
resorting to the asymptotic approximations. Moreover, 
we think that, in applications, continuous-time random 
walks are seen as a more fundamental model than frac- 
tional diffusion equations, and sample paths will be gen- 
erated without taking the scale parameters 7^ and 74 to 
the diffusive limit, by using the approach presented in 
this paper. 
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